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SYMBOLS 


Number  of  input  samples 
Number  of  partial  transformations 
Continuous  -  frequency  variable 
Highest  frequency  in  a  spectrum 
Sampling  rate  or  frequency 

Frequency  increment  between  successive  components 
Continuous  time  variable 

Effective  period  for  a  time  function  when  it  is  periodic 
Time  increment  between  successive  samples 
The  set  of  input  data  samples 

Where  £-0,  1,  2,  ....  N-l  and  m-1,  2,  . . .M  is  an  output 
sequence  of  the  mth  nodal  column  computations 

Becomes  the  input  array  during  the  (nr4-l)st  stage  of  computations 
where  Xjj+iCi)  is  the  output  array 

Fundamental  frequency 

Number  of  samples  per  cycle  of  f(t) 


I.  INTRODUCTION 


In  recent  years  the  use  of  digital  simulation  has  become  an  intricate 
part  of  military  and  civilian  projects.  A  significant  number  of  these  projects 
involve  digital  system  simulations  which  require  frequency  response  analysis 
of  the  ensuing  discrete-time,  discrete-frequency  simulation  output.  Very 
often  these  outputs  are  not  composed  of  simple  frequency  sinewaves,  but  rather 
contain  extraneous  harmonics  along  with  the  desired  fundamental  frequency, o>0  » 
making  it  impractical  to  calculate  amplitude  ratios  over  the  desired  frequency 
range.  When  this  is  the  case,  the  frequency  response  of  the  system  can  be 
obtained  from  the  input-output  amplitude  ratios  of  the  input  and  output 
frequency  spectrums  at  each  fundamental  frequency,  o)0  ,  over  the  desired 
frequency  range. 

An  ideal  method  of  generating  the  frequency  spectrum  of  a  discrete-time, 
discrete-frequency  signal  of  a  computer  simulation  is  a  Fast  Fourier  Transform 
(FFT)  on  the  simulation  output  data.  The  FFT  generates  a  discrete  frequency 
spectrum  analogous  to  the  Fourier  spectrum.  Both  give  two  impulses  at  +  o>o 
and  - <d0  for  a  sinewave  input.  The  frequency  response  of  a  discrete  system 
simulation  can  be  found  by  inputting  a  sinewave  into  the  discrete  system 
simulation  and  using  the  frequency  spectrum  iftput  and  output  Impulse  ratios 
at  oj0  over  the  desired  frequency  range. 

A  FFT  computer  subroutine  using  VAX-11  FORTRAN  has  been  written  to  perform 
the  FFT.  The  subroutine  can  be  linked  with  a  system  simulation  to  provide  the 
frequency  spectrum  Impulse  data  as  a  part  of  the  system  simulation.  This 
report  outlines  the  development  and  checkout  of  the  FFT  routine. 

II.  BACKGROUND 

The  FFT  is  an  algorithm  for  computing  the  discrete  Fourier  transform  of 
discrete  data  samples.  There  are  many  available  FFT  algorithms.  The  Radix-2 
FFT  is  the  one  most  commonly  used.  It  is  based  on  representing  an  array  of 
size  N“2M  as  a  product  of  M  factors,  each  of  which  is  equal  to  2.  Radix-2 
FFT  algorithms  are  derived  by  decomposing  the  discrete  Fourier  transform  into 
successively  smaller  discrete  Fourier  transforms.  The  manner  of  the 
decomposition  produces  the  variation  found  in  Radix-2  algorithms.  Most  of 
these  algorithms  may  be  classified  as  follows: 

A.  Decimation  in  Frequency 

1.  In-place  algorithm  or 

2.  Natural  input-output  algorithm 

B.  Decimation  in  Time 

1.  In-place  algorithm  or 

2.  Natural  input-output  algorithm 


The  FFT  algorithm  used  to  write  the  FFT  subroutine  is  an  in-place 
decimation  in  frequency,  Radix-2  algorithm  originally  proposed  by  Gentlemen 
and  Sande. 


The  algorithm  is  implemented  in  the  following  steps: 

a.  Initialization 

b.  A  sequence  of  transformations,  one  partial  transformation  for 


each  factor. 

c.  An  unscrambling  procedure. 

III.  DERIVATION  OF  A  RADIX-2  FFT  BY  DECIMATION  IN  FREQUENCY 

The  discrete  Fourier  transform  of  {  x(n)  |  is  a  periodic  sequence  of  complex 
numbers  |  X(k),k  -  0,1,...,  N-l  ^  ,  defined  by 


X(k) 


N-l 
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n-0 


x(n)  W 


nk 


(1) 


where 


W  -  exp  (-J-tt)  -  cos  (^-)  -j  sin  (~) 


Dividing  the  input  sequence  into  two  halves  gives 
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Combining  the  two  summations  in  equation  3  and  using  the  fact  that 

(!)k  /  xk 

2  -  (-l)  ,  yields 


(3) 


X(k) 


(x(n)  +  (-1)  x  (n  +  J  w 


(4) 


n**o 

Since  (-l)^  is  equal  to  1  for  even  K  and  equal  to  -1  for  odd  K,  let 
even  K  *  2r  and  odd  K  ■  2r  »  1. 

Dividing  the  sequence  in  equation  4  into  an  even  and  odd  sequence  yields 
a  decimation  in  frequency  of 

(f)-l 

X(2r)  *  ^  *  ^x(n)  +  x(n  +  7  ))  W  2rn 

n=o 

(f  )  -1 

X(2r  +1)  - 

n*o 

The  arrays  X(2r)  and  X(2r  +  1)  are  the  7  points  DFT  of  the  input  arrays 
shown  in  equation  (5).  The  signal  flow  graph  for  an  eight-point  input  (N=8) 
is  shown  in  Figure  1. 


^x(n)  - 


x(n  + 


)^j  W°  w 2 


an  N-point  DFT  computation  into  two  N/2-point  DFT  computations 


The  bit-reyejcsed  output  can  also  be  determined  by  the  following  method; 


A,  Write  out  the  sequence  of  integer  index  numbers, 

B,  Convert  the  sequence  to  binary  form. 

C,  Reverse  the  order  of  all  the  bits  in  each  of  the  binary  index 

numbers . 


D.  Convert  back  to  integer  index  numbers. 

Table  1  shows  the  determination  of  the  bit-reversed  output  for  the  N=8 
example  using  the  method  just  described. 


TABLE  1,  Index  Integers  and  Their  Bit-Reversed  Output 
Integers  For  N«8 


■>,!.rn>.iTmj 
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6 
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Binary  Index 

000 

001 

010 

011 

100 

101 

110 

111 

Binary 

Bit -Reversed 

000 

100 

010 

110 

001 

101 

011 

111 

Bit-Reversed 
Index  Integer 

H 

m 

2 

6 

1 

5 

■1 

7 
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Repetition  of  this  decomposition  for  N=8  leads  to  Figures  2  and  3. 


Figure  2. 


Figure  3. 


Flow  graph  of  decimation-in-frequency  decomposition  of  an 
eight-point  DFT  into  four  two-point  DFT  computations . 


Flow  graph  of  complete  decimation- in- frequency  decomposition 
of  an  eight-point  DFT  computation. 


The  preceding  Figures  show  that  the  input  sequence  is  in  a  natural  order 
and  the  output  sequence  is  in  an  unnatural  order.  This  input-output  order  is 
referred  to  as  natural  input,  bit-reversed  output.  It  is  a  result  of  decom¬ 
position,  Recall  that  originally  the  input  sequence  is  divided  into  even- 
numbered  samples  and  odd-numbered  samples  with  the  even-numbered  samples  in 
the  first  half  and  the  odd-numbered  samples  in  the  second  half.  Such 
separations  are  carried  out  by  examining  the  lease  significant  bit  of  the 
binary  index  representation  where  a  zero  corresponds  to  an  even  number  and 
a  one  corresponds  to  an  odd  number,  When  the  even  and  odd  subsequences 
are  sorted  into  their  even  and  odd  parts,  the  second  least  significant  bit  of 
the  binary  index  representation  is  examined.  This  process  is  repeated  until 
N  subsequences  of  length  1  are  obtained.  This  sorting  process  for  N  =  8  is 
shown  in  Figure  4. 


x(000)  0 

xOOO)  4 
x(OlO)  2 
*(110)  6 
*(001)  1 
*(101)  5 

x(0 1  i )  3 
*(i  1  l)  7 


Figure  4.  Bit-reversed  sorting  for  N=8. 


IV.  COMPUTATIONS 


The  vertical  nodes  in  the  preceding  flow  graphs  correspond  to  successive 
inner  column  computations.  Each  inner  column  computation  takes  a  set  of  N 
complex  numbers  and  transforms  them  into  another  set  of  N  complex  numbers. 
This  process  is  repeated  M*  log2  N  times. 

The  flow  graph  for  the  basic  computation  (the  butterfly  computation)  is 
shown  in  Figure  5. 


The  set  of  equations  found  from  the  flow  graph  shown  in  Figure  5  is  as 
follows 


Vi  <»>  '  x„  +  x„  <«> 
Vi  (<1)  ■  x*  <»>  -  x»  M 


where 


tn  ■  0,  1,  . . .  M  -  1 

p  -  1,  2,  ...  N/2 
N 

q  "  2  +  1,  ...  N 


These  equations  are  the  algorithm  for  the  transformations  necessary  for 
the  Radix- 2-  FFT.  Looking  at  equation  (6),  the  data  pair  (^(p),  X^q))  is 
used  once -to  compute  the  new  data  pair  (XnrH  (p) ,  3Cnrfi(q))  and  is  not  used 
again.  Thus  the  transformations  defined  by  equation  (6)  can  be  performed  in 
place  where -each  new  transformation  is  stored  in  the  same  location  occupied 
by  the  preceding  transformation.  This  can  be  seen  graphically  in  Figure  3 
by  the  horizontal  lines  which  connect  consecutive  nodes. 
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•V 


As  an  example  of  the  use  of  the  set  of  transformation  equations,  the 
first  Inner  column  computations  of  the  8-point  FFT  as  shown  in  Figure  5 
are  found  to  be: 


Xj  (0) 
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X  (0)  +  X  (4) 

o  o 

Xl  (1) 

- 

X  (1)  +  X  C5) 
o  0 

Xx  (2) 

» 

X  (2)  +  X  (6) 

o  0 
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- 
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o  0 
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m 
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Xx  W 

=* 

w2  [xo(2)  -  X0(6)] 
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m 
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The  other  column  computations  can  be  found  in  a  similar  manner  using 
the  transformation  algorithm  as  shown  in  equation  (6) .  The  number  of  column 
computations  will  always  be  equal  to  M. 


V.  PARAMETERS 


Selection  of  FFT  Parameters: 

A.  The  number  of  input  samples  (N)  must  be  a  power  to  two. 

B.  The  sampling  rate  (fs)  must  be  greater  than  twice  the  highest  possible 
frequency  in  the  spectrum  (fh)  to  avoid  aliasing. 

C.  The  period  of  the  time  function  (tp)  is  chosen  to  be  longer  than  the 
time  length  of  the  signal  to  prevent  overlapping. 

D.  The  desired  frequency  resolution  (F)  is 


E.  The  only  way  to  select  both  tp  and  fs  independently  is  to  use  the 
relationship. 


N  > 


(8) 


VI.  RADIX- 2  FFT  COMPUTER  PROGRAM 

A  computer  program  which  uses  the  Radix-2  FFT  algorithm  developed  in  this 
report  and  containing  a  binary  unscrambling  algorithm  is  shown  in  the  Appendix. 

The  computer  program  is  written  as  a  subroutine.  There  are  three  inputs 
necessary  to  run  the  FFT  subroutine. 

A.  N 

B.  M 

C.  The  input  samples. 

The  input  samples  (N)  are  placed  into  two  arrays  as  listed  below. 

P  (N,  1)  Real  part  of  sample 

P  (N,  2)  Imaginary  part  of  sample 

If  the  input  samples  are  all  real  numbers,  F  (N,2)  would  be  zero 

filled. 


In  order  to  Insure  reliable  FFT  results,  care  should  be  taken  in  selecting 
the  input  data  sampling  rate  and  sampling  duration. 

The  output  of  the  FFT  subroutine  is  contained  in  the  two  arrays  P  (N,  1) 
and  P  (N,  2) .  The  series  of  magnitudes  found  from  these  arrays  is  the 
discrete  Fourier  series  equally  spaced  F  units  apart.  Recalling  that  a 
Fourier  series  separates  a  periodic  function  of  period  T  into  sinuisoldal 
components  of  frequency  (o0  ,  2<jj0  ....  .r^o  *  where  0)o*2tt/t  *8  the  fundamental 
frequency  and  the  other  frequencies,  2oJq  , . . .  ,rw0 ,  are  the  harmonics  of  <o0, 
the  output  of  the  FFT  subroutine  is  in  fact  a  harmonic  analysis  where  the 
output  magnitudes  are  amplitudes  of  signal  components  are  discrete  frequency 
intervals . 

VII.  CHECKOUT 

The  FFT  program  was  checked  out  by  inputting  samples  from  the  function 


f(t)  ■  5  sin  <Ujt  +  10  sin  co2t 


1th 


6.28  RAD/ SEC 

0>2  -  31.42  RAD/SEC 

1.0  Hz  and 

f2  -  5  Hz 

The  following  FFT  parameters  were  used: 

N  -  1024 

tn  * 


2.0  sec 


The  output  of  the  FFT  program  is  shown  in  Table  2.  The  table  shows 
the  expected  results  of  a  5  degree  signal  at  1  Hz,  and  a  10  degree  signal 
at  5  Hz.  -  which  are  the  same  results  that  would  be  obtained  from  a  DFT. 


TABLE  2.  FFT  Output  For  The  Checkout  Example 


FREQ 

SPECTRUM  MAGNITUDE 

(Hz) 

(DEGREE) 

■j  . 

3.3  377379E-35 

^  t  3  v  'J  «  C  w  C 

7.535c ‘♦4  9E-.4 

i. or  ococ 

s.'w  3C  381 

1 • 50 vC >  0 

9.  o5  3E  331  £-:. 4 

2.  10. Cj u 

1  •  -»3  37  9 2 7E -C  3 

C  .  S0-!'«  jw 

2.4941273E-. 

3.0C0vji 

1.7234  333:-.  3 

3  .  5  Cj  C  j  ’ 

1.2C 0;881crv 3 

*• J 0  j u 

1.9S7S  39  3E-.  3 

h • 50  jOw  w 

1.1196  ’62E-C  3 

3  .  J  v  C  .  U 

9.999151 

*  5  •  5  *  w  c  0  » 

9. 76  3<.  62  7E-*.  4 

5.1663137E-.4 

6  •  5w00Jj 

1 .  390C  615E-. 3 

7.  j  C  ;G  7  0 

3.2764  J91B-v4 

7  •  5  C  V  C  *  0 

9.  )2  8E431E-C.  4 

d  •  v  t>  0  C  • '  0 

7. 7732962c-. 4 

3  •  5  C  0  *  0  « 

7.4152361c-'.  4 

9 . J  0  J  0 .0 

1 .  5  92  474  c-'.  3 

9a  SC  vC 10 

3.626S  751E-C4 

1 «/  •  0  }  C  J  j 

1.71165 v 6c- 13 

To  check  the  effect  of  sampling  rate  and  the  length  of  record  (tp). 
Table  3  was  constructed.  A  spectral  analysis  (FFT)  of  f(t)  *  10  cos  20  irt 
was  done  for  the  cases  tp  ■  .1,  1,  and  4  sec.  and  f8  was  varied  until  the 
minimum  sampling  frequency  was  found. 

Table  3  shows  that  a  minimum  sampling  frequency  for  the  10.0  Hz  signal 
is  256  Hz  for  tp  equal  to  1  and  4.  The  sampling  rate  at  tp  equal  to  .1  is 
larger,  but  this  is  due  to  the  discrete  values  that  N  must  have  which  then 
causes  a  discrete  jump  in  the  sampling  frequency  as  seen  between  tp  »  .1  and 
1  sec.  The  minimum  sampling  frequency  produces  an  impulse  at  oj0  and  zeroes 
at  all  other  frequencies.  It  should  be  noted  that  at  lower  sampling  frequen¬ 
cies  the  impulse  at  oj0  can  be  seen  along  with  other  low  impulses  at  other 
frequencies. 


TABLES  3.  Minimum  Sampling  Times  and  Frequency  for 
f(t)  ■  10  cos  20  irt . 


.1 

1 

640 

256 

64 

256 

1 

10 

The  main  sources  of  error  for  the  FFT  are  sampling  rate,  quantization 
and  round-off.  The  sampling  rate  error  as  given  by  Hopper  and  Newberry  [3.] 
is 


E  <  — - - 

s  —  12  K2  percent 


where  K  is  the  number  of  samples  taken  ®n  each  cycle  of  f(t),  and  the 
quantization  and  round-off  error  as  given  by  Welch  [4.]  is 


Eq  <  0.1  percent  (10) 


Using  equation  (9)  with  the  data  in  Table  3.,  Es  is  found  to  be  zero  for  t 

equal  to  .1  sec, and  0.01  percent  for  t  equal  to  1  and  4  sec.  r 

•  P 

VIII.  CONCLUSION 

The  FFT  subroutine  presented  in  this  report  gives  good  results  in  the 
frequency  domain.  The  results  are  discrete  and  are  nonexistant  between 
adjacent  frequency  intervals. 

Nothing  is  known  between  the  discrete  frequencies.  For  more  resolution 
of  frequency,  the  record  lenght,  tp,  must  be  increased.  Care  must  be  taken 
to  ensure  that  co0  is  one  of  the  discrete  frequencies  calculated,  and  the 
sampling  frequency  is  high  enough  to  yield  low  error  results. 

Overall,  the  FFT  subroutine  is  a  very  useful,  fast  computational 
algorithm  which  can  be  used  with  any  digital  system  simulation  when  frequency 
spectrum  processing  is  needed  in  the  calculation  of  the  system's  frequency 
response . 
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SUBROUTINE  FFT 
COMMON  C (2000) 

EQUIVALENCE  (C(  900), IP9  ) 
EQUIVALENCE  (C(  901), IP8  ) 
DIMENSION  P (8 192, 2) 


PARAMETERS 

INPUT***********************FOURIER  TRANSFORM************************ 

P  -  AN  ARRAY  DIMENSIONED  MAX — P  (8192,2) 

DIMENSIONED--? (IP 9, 2)  USED  FOR  THE  DATA 
IP9  -  N  «  A  POWER  OF  2  -  THE  NUMBER  OF  DATA  POINTS 
IP9  MUST  BE  A  POWER  OF  2 

IP9  *  2**IP8  IP8  =•  M  MAX - 8192  -  2**13 

OUTPUT**********************FOURIER  TRANSFORM************************ 

P  -  THE  TRANSFORMED  DATA.  THIS  MEANS  THE  INITIAL  DATA  IS  DESTROYED  BY 
THE  ROUTINE  AND  REPLACED  WITH  THE  TRANSFORMED  DATA. 


*****ppT***ft* 

IQ1  -  2*IP9 
Q6-3.141592654/IP9 
DO  30  T-3r,IP8 
IQ1-IQ1/2 
IQ2-IQ1/2. 

Q6-Q6*2.0 

JQ“IP9+1-IQ1 

DO  20  J-1,JQ,IQ1 

Q5— Q6 

IP4-J-1 

DO  10  K-1,IQ2 

Q5-Q5+Q6 

IQ3-IP4-HC 

IQ4-IQ3+IQ2 

Q7-COS(Q5) 

Q8-SIN(Q5) 

Q9~Q7* (P (IQ3 , 1 )-P (IQ4 , 1 ) )+Q8* (P (IQ4 , 2) -P (IQ3 , 2) ) 
Q8-Q7* (P (IQ3, 2)-P (IQ4 , 2) )+Q8* (P (IQ3 , 1 ) -P (IQ4 , 1) ) 
P (IQ3 , 1)*P (IQ3 , 1)+P (IQ4 , 1 ) 
P(IQ3,2)-P(IQ3,2)+PIQ4,2) 

P(IQ4,1)-Q9 
P(IQ4,2)-Q8 
10  CONTINUE 
20  CONTINUE 
30  CONTINUE 

DO  60  1-1 ,  IP9 

IQ3-I-1 

IQ4-0.0 


A-l 


DO  50  J-1.IP8 
IQ4-2*IQ4 
P4-FL0AT(IQ3)/2,0 
IF (INT(P4) ,EQ,P4)GOTO  40 
IQ4-IQ4+1 
40  IQ3-INT(P4) 

50  CONTINUE 


IQ4-IQ4+1 

IF(IQ4.LE.I)G0T0  60 
Q9-P (I , 1 ) 

Q8-P(I,2) 
P(I,1)«P(IQ4,1) 
P(I,2)-P(IQ4,2) 
P(IQ4,1)-Q9 
P(IQ4,2)-Q8 
60  CONTINUE 


DO  70  I-1.IP9 
P(I*1)“P(I»1)/IP9 
P(I,2)-P(I,2)/IP9 
70  CONTINUE 


J-IP9 

K-IP9/2 

DO  80  1-2, K 

QI-P(I,1) 

Q2-P(I,2) 

P(I,1)-P(J,1) 

P(I»2)*P (J,2) 

P(J,1)-Q1 

P(J,2)-Q2 

J-J-l 

80  CONTINUE 


PRINT  OUT  FREQ  (FQFFT)  AND  MAG  (FFTM)  OF  FFT  *************** 
DO  90  1-1,30 

FF1W-SQRT ( (P  (1 , 1 ) *2 . 0) **2+(P ( 1 , 2 ) *2 . 0) **2) 

PRINT*, FQFFT, FFTM 
FQFFT-FQFFT+1 . 0/ (IP9*DTFFT) 

90  CONTINUE 

END  IF 

RETURN 

21  FORMAT(2(lx,G14.6)) 
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